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Donor— acceptor  binding  of  the  jr-electron-poor  cyclophane  cyclobis(paraquat-/?-phenylene)  (CBPQT4+)  with 
the  jr-electron-rich  tetrathiafulvalene  (TTF)  and  1,5-dioxynaphthalene  (DNP)  stations  provides  the  basis  for 
electrochemically  switchable,  bistable  [2]rotaxanes,  which  have  been  incorporated  and  operated  within  solid- 
state  devices  to  form  ultradense  memory  circuits  ( ChemPhysChem  2002,  3,  519—525;  Nature  2007,  445 , 
414—417)  and  nanoelectromechanical  systems.  The  rate  of  CBPQT4+  shuttling  at  each  oxidation  state  of  the 
[2]rotaxane  dictates  critical  write-and-retention  time  parameters  within  the  devices,  which  can  be  tuned  through 
chemical  synthesis.  To  validate  how  well  computational  chemistry  methods  can  estimate  these  rates  for  use 
in  designing  new  devices,  we  used  molecular  dynamics  simulations  to  calculate  the  free  energy  barrier  for 
the  shuttling  of  the  CBPQT4+  ring  between  the  TTF  and  the  DNP.  The  approach  used  here  was  to  calculate 
the  potential  of  mean  force  along  the  switching  pathway,  from  which  we  calculated  free  energy  barriers. 
These  calculations  find  a  turn-on  time  after  the  rotaxane  is  doubly  oxidized  of  ~10-7  s  (suggesting  that  the 
much  longer  experimental  turn-on  time  is  determined  by  the  time  scale  of  oxidization).  The  return  barrier 
from  the  DNP  to  the  TTF  leads  to  a  predicted  lifetime  of  2.1  s,  which  is  compatible  with  experiments. 


1.  Introduction 

The  electrochemically  switchable,  bistable  [2]rotaxanes1 
(Figure  1)  developed  in  recent  years  by  Stoddart  and  co-workers 
exhibit  two  distinct  co-conformations:3-7  the  ground-state  co¬ 
conformation,  in  which  the  cyclobis(paraquat-p-phenylene) 
(CBPQT4+)  encircles  the  tetrathiafulvalene  (TTF)  station,  and 
the  metastable  state  co-conformation,  in  which  the  CBPQT4+ 
encircles  the  1,5-dioxynaphthalene  (DNP)  station.2,8-15  The 
population  of  the  two  co-conformations  may  be  shifted  away 
from  equilibrium  by  temporarily  oxidizing  one  or  two  electrons 
from  the  TTF  units.  This  switching  process  forms  the  basis  of 
using  these  compounds  as  storage  elements  in  molecular 
electronic  devices.  Consequently,  significant  experimental  efforts 
have  been  made  to  investigate  the  switching  behavior  of 
molecular  switches5,6,16-23  and  molecular  machines24-29  in 
various  environments,  such  as  solution,5,6,30-37  polymer  elec¬ 
trolyte  gels,38  metal  surfaces39-41  and  devices.7,16,17  Important 
experimental  evidence15  for  molecular  switching  in  these  devices 
was  the  correlation  of  the  kinetics  of  relaxation  from  the  DNP 
to  the  TTF,  across  each  of  these  environments.  However,  the 
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rate  of  this  process  is  also  a  function  of  the  molecular  structure, 
suggesting  that  longer  storage  times,  or  even  nonvolatile 
memory,  might  be  possible  with  the  appropriate  molecular 
design. 

Computation  chemistry  calculations  could  provide  an  effective 
approach  for  optimizing  the  performance  of  such  molecular 
switches,  but  such  applications  require  that  the  accuracy  of  the 
theory  be  validated  by  comparing  to  well-documented  experi¬ 
mental  results.  The  purpose  of  this  paper  is  to  provide  such 
validation.  Previously,  we  investigated  these  compounds  using 
a  multiscale  first  principles  approach  combining  quantum 
mechanics  (QM)  and  atomistic  force  field  (FF)  methods  42-47 
First  we  considered  the  molecules  as  individual  species,  and 
then  we  examined  self-assembled  monolayers  bound  to  gold 
surfaces  or  compressed  into  Langmuir  monolayers  at  the 
air— water  interface.  These  studies  successfully  predicted  a 
number  of  phenomena  that  were  confirmed  later  experimentally, 
including  the  higher  conductivity47  of  the  DNP  relative  to  the 
TTF,  and  the  increased  stability  of  the  TTF  relative  to  the  DNP 
(by  2.0  kcal/mol  from  QM,  2.3  kcal/mol  from  the  FF,  and 
1.4— 1.6  kcal/mol  from  experiment).15,43  In  addition,  on  the  basis 
of  the  predicted  footprint  of  the  115  A2/molecule  for  the  self- 
assembled  structure,  we  predicted  that  the  surface  tension  of 
the  TTF  is  32%  lower  than  that  of  the  DNP,  an  observation 
that  was  confirmed  in  subsequent  experiements  43,44 

In  this  study,  we  evaluated  the  free  energy  profile  of  the 
shuttling  motion  of  the  CBPQT4+  ring  between  the  TTF  and 
the  DNP  stations  to  determine  how  the  nature  of  the  rotaxane 
affects  the  switching  and  relaxation  rates.  These  rates  have  been 
determined  experimentally  in  various  environments,15,35,40,48  and 
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Ground  State  Metastable 

TTF  DNP 

Figure  1.  Structural  formula  of  the  two  co-conformations  of  a  bistable  [2]rotaxane  fragment  used  in  this  study. 


we  now  want  to  understand  the  structural  contributions  to  these 
rates.  We  seek  to  find  a  level  for  theoretical  calculations  of  these 
rates  that  is  both  accurate  and  fast  so  that  we  can  use  theory  to 
optimize  the  structural  characteristics  to  achieve  desired  rates. 

Rather  than  finding  the  energy  barrier  for  a  minimized 
reaction  path  connecting  the  two  states,  we  used  potential  of 
mean  force  (PMF)  to  evaluate  the  change  of  free  energy  along 
the  shuttling  pathway  of  the  CBPQT4+  ring  between  the  TTF 
and  the  DNP  so  that  we  can  determine  the  rates  at  the 
experimental  temperature.  We  carried  out  these  calculations  for 
three  oxidation  states  of  the  molecule  relevant  to  the  switching 
and  thermally  activated  relaxation  process. 

2.  Simulation  Details 

2.1.  Potential  of  Mean  Force  from  Constrained  Molecular 
Dynamics  Simulation.  The  experimental  time  scale  for  the  ring 
to  relax  back  from  the  DNP  to  the  TTF  is  lO’-lO3  seconds,1534,49 
suggesting  that  simple  molecular  dynamics  (MD)  simulation 
of  a  few  hundred  nanoseconds  might  not  be  sufficiently  ergodic 
to  provide  an  accurate  transition  rate.  Hence,  we  adopted  the 
“Blue  Moon  sampling”  technique50,51  of  constrained  MD 
simulations  using  holonomic  constraints  that  fix  the  systems 
along  the  reaction  coordinate.  To  determine  the  free  energy 
barrier,  we  used  the  reaction-coordinate  (/^-dependent  potential 


of  mean  force  (PMF),  Frxn(R)  defined  as  the  integration  of  the 
mean  force  (MF)  along  the  reaction  coordinate,  —  dFrxn(R)/dR.52 

pR  d Frxn(R') 

=  Fm(oo)  +  l  dR'  (1) 

Here,  the  MF  is  a  measurable  quantity  from  our  simulations. 
To  calculate  the  MF,  we  assumed  that  the  CBPQT4+  ring  moves 
between  the  TTF  and  the  DNP  along  the  backbone  of  the 
rotaxane  (Figure  2a),  which  we  assume  to  be  in  an  extended 
conformation  but  with  the  minimized  structure.  This  extended 
conformation  should  provide  the  fastest  shuttling  motion  of  the 
CBPQT4+  ring,  being  governed  mainly  by  its  interaction  with 
the  backbone.  This  MF  does  not  account  for  the  presence  of 
folded  chain  conformations,  so  that  the  PMF  may  lack  some 
contributions  from  conformational  entropy. 

First,  we  prepared  the  extended  rotaxane  backbone  without 
the  CBPQT4+  ring  using  quantum  mechanical  geometry  opti¬ 
mization  at  the  level  of  B3LYP/6-31G*  (Figure  2a).  Then,  we 
added  and  optimized  the  CBPQT4+  ring  at  various  fixed  points 
on  the  fixed  extended  backbone  (Figure  2b)  using  quantum 
mechanics.  Thus,  the  atomic  partial  charges  of  ah  atoms  are 
allowed  to  readjust,  depending  on  the  relative  position  of  the 
charge  acceptor  (CBPQT4+)  with  respect  to  the  charge  donor 
(TTF  and  DNP). 
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Figure  2.  (a)  Backbone  of  the  rotaxane  molecule  simulated  in  this 
study,  (b)  CBPQT4+  ring  positions  along  the  backbone  (unit:  A). 


Figure  3.  Charges  for  the  initial  nine  structures  obtained  from  QM 
with  Mulliken  analysis.  In  addition,  we  included  two  more  structures 
beyond  each  station  of  the  TTF  and  the  DNP,  using  charges  identical 
to  those  for  the  equilibrium  CBPQT4+@TTF  and  CBPQT4+@DNP 
cases,  respectively.  Ten  more  structures  were  generated  on  the  basis 
of  these  eleven  structures.  The  position  of  the  CBPQT4+  ring  for  each 
additional  structure  was  obtained  using  the  arithmetic  average  of  the 
two  adjacent  cases  from  the  eleven  structures.  The  charges  were  also 
averaged. 


To  obtain  the  change  of  the  PMF  during  the  shuttling  process, 
we  first  evaluated  the  MF  as  a  function  of  the  position  of  the 
CBPQT4+  ring.  Because  the  length  of  the  backbone  is  54.7  A 
and  the  distance  between  the  TTF  and  the  DNP  is  36.9  A,  we 
chose  to  sample  the  dynamics  for  nine  independent  samples, 
each  of  which  has  the  z-coordinate  (along  the  backbone)  of  the 
center  of  mass  (COM)  of  the  CBPQT4+  ring  at  a  different 
position  along  the  extended  backbone,  as  schematically  pre¬ 
sented  in  Figure  3.  Using  quantum  mechanics,  the  geometry 
and  atomic  charges  were  obtained  from  each  of  these  nine  cases. 

After  preparing  these  nine  initial  structures,  we  prepared  two 
more  structures  beyond  each  station  of  the  TTF  and  the  DNP 
with  identical  charges  to  the  CBPQT4+@TTF  case  and  the 
CBPQT4+@DNP  case,  respectively.  In  addition,  we  constructed 
another  ten  structures  in  which  the  position  and  charges  of  the 
CBPQT4+  ring  were  calculated  by  arithmetically  averaging  the 
coordinates  and  charges  of  two  consecutive  structures  in  the 


eleven  structures.  Thus,  a  total  of  21  structures  were  prepared 
for  simulations. 

Then,  to  simulate  both  the  turning  on  and  turning  off  the 
rotaxane  switch,  we  investigated  the  effect  of  oxidation  of 
rotaxane  molecule  on  the  free  energy  profile,  for  three  different 
oxidation  states:  the  neutral  state,  the  +1  oxidation  state,  and 
the  +2  oxidation  state. 

The  QM  calculations  of  the  charges  for  the  nine  different 
structures  were  repeated  for  each  of  the  three  oxidation  states: 
0,  +1,  and  +2.  The  atomic  partial  charge  distributions  are 
tabulated  in  the  Supporting  Information,  Table  SI— S3. 

All  quantum  mechanical  computations  in  this  study  were 
performed  using  Jaguar.53 

2.2.  Constrained  Molecular  Dynamics  Simulation.  Next, 
we  carried  out  a  constrained  NVT  MD  simulation  at  300  K  for 
500  ps  to  equilibrate  each  system.  This  MD  was  then  continued 
for  an  additional  3  ns  at  300  K  (constrained  NVT  MD)  to 
compute  the  MF.  The  constraint  was  introduced  using  Gauss’ 
principle  of  least  constraints54  to  fix  only  the  z-component  of 
the  center  of  mass  (COM)  of  the  CBPQT4+  ring  parallel  to  the 
molecular  axis  direction  (z-axis  direction  as  in  Figure  2b).  To 
ensure  that  our  constraint  dynamics  produces  the  correct 
equilibrium  averages  without  bias  due  to  ensemble  sampling, 
we  used  Fixman’s  theorem55  to  evaluate  the  metric  effect 
originating  from  the  holonomic  constraints.  We  determined  that 
the  metric  effect  only  adds  a  constant  scalar  value  to  the  absolute 
free  energy  values,  which  has  no  influence  on  the  relative 
energetics.  (Details  are  in  the  Supporting  Information.) 

We  also  fixed  the  position  of  the  last  oxygen  atom  at  each 
end  of  the  backbone  to  retain  the  extended  conformation.  This 
restricts  the  conformational  flexibility  of  the  system,  which 
suppresses  conformational  entropic  contributions  to  the  free 
energy.  The  mean  force  was  sampled  from  such  constrained 
MD  simulations. 

2.3.  Force  Field  and  MD  Parameters.  We  used  the  generic 
DREIDING  force  field,56  which  was  found  to  lead  to  accurate 
results  in  our  previous  studies  on  rotaxane  systems  43-45  It  was 
also  successful  in  our  studies  on  various  other  molecular 
systems,  such  as  the  hydrated  polymer  electrolyte  membranes57-59 
and  the  surfactant- mediated  air— water  interface.60,61 

The  force  field  has  the  form 

^total  ^vdW  ^bond  ^angle  ^torsion  ^inversion  (^) 

where  Ftotai?  ^vdw?  Eq?  Ebond?  Eang\e,  Etorsion?  and  Ei nversion  &re  the 
total  energies,  the  van  der  Waals,  electrostatic,  bond  stretching, 
angle  bending,  torsion,  and  inversion  energy  components, 
respectively,  and  the  force  field  parameters  are  described  in  the 
original  papers.56  The  atomic  charges  were  obtained  from  a  QM 
Mulliken  population  analysis  as  indicated  above. 

All  MD  simulations  were  performed  using  LAMMPS  (large- 
scale  atomic/molecular  massively  parallel  simulator)  MD  code 
from  Plimpton  at  Sandia. 62,63  The  equations  of  motion  were 
integrated  using  the  velocity— Verlet  algorithm,64  with  a  time 
step  of  0.01  fs.  This  unusually  small  time  step  was  to  ensure 
high  quality  sampling  of  phase  space  by  avoiding  abrupt  changes 
in  atomic  positions. 

The  temperature  was  kept  constant  during  the  MD  using  the 
Berendsen  thermostat  with  temperature  damping  time  of  0.01 
fs.  To  demonstrate  that  our  MD  leads  to  a  proper  canonical 
ensemble,  the  probability  distribution  function  (PDF)  of  kinetic 
energy  KE  (=mv1/ 2)  is  shown  in  Figure  4.  The  PDF  is  quite 
close  to  the  Maxwell— Boltzmann  distribution  of  energy  at  T  = 
300  K,  indicating  that  the  simulation  describes  a  proper 
canonical  ensemble.  Furthermore,  the  PDF  for  each  component 
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Figure  4.  Probability  density  function  of  kinetic  energy  KE  (=mv2/ 2) 
is  from  the  MD  simulation  of  the  CBPQT4+  ring@TTF  (blue  line)  at 
300  K.  Here  the  time  step  was  0.01  fs  and  the  total  simulation  time 
was  3  ns  after  500  ps  of  equilibration.  The  black  dashed  line  compares 
with  the  Maxwell— Boltzmann  distribution  of  the  energy,  2(KE/ 
jt(kBT?)m  Qxp(-KE/kBT),  for  T  =  300  K. 


Simulation  time  (ps) 

Figure  5.  Change  of  mean  force  as  a  function  of  simulation  time.  In 
each  case  this  follows  500  ps  of  equilibration  time.  This  plot  shows 
two  representative  cases:  the  CBPQT4+  ring@TTF  (TTF)  and  the 
CBPQT4+  ring@DNP  (DNP). 


of  velocity  is  the  same  and  the  system  obeys  the  equipartition 
theorem  (see  Figure  SI  in  the  Supporting  Information). 

Figure  5  shows  the  typical  behavior  of  the  MF  as  a  function 
of  simulation  time  for  two  representative  systems:  one  is  the 
ground  state,  CBPQT4+@TTF,  green  color,  denoted  as  TTF  and 
the  other  is  the  metastable  state,  the  CBPQT4+  ring  on  the  DNP, 
red  color,  denoted  as  DNP.  This  shows  that  the  mean  force  was 
well  equilibrated  for  both  cases. 

The  weakness  of  this  blue  moon  sampling  method  is  that  the 
error  in  each  MF  measurement  is  integrated  to  obtain  the  PMF 
profile  along  the  reaction  coordinate.  From  block  averages,  we 
estimate  the  uncertainty  of  the  MF  values  to  be  0.04  (kcal / 
mol)/A  for  CBPQT4+@TTF  and  0.22  (kcal/mol)/A  for 
CBPQT4+@DNP.  Assuming  that  these  errors  are  random  and 
that  the  average  value  is  0.13  (kcal/mol)/A,  we  estimate  that 
the  error  of  the  free  energy  difference  between  two  stations  is 


Figure  6.  Change  of  potential  of  mean  force  as  a  function  of  ring 
position  along  the  backbone.  The  blue  curve  allows  the  charge  to  change 
adiabatically  as  the  ring  moves  along  the  dumbbell,  which  is  the  reliable 
result.  The  other  two  curves  show  the  error  obtained  when  the  charges 
are  fixed:  the  green  curve  uses  fixed  charges  from  the  ring  @ TTF;  the 
red  curve  uses  a  fixed  charge  from  the  ring  @  DNP. 


38.8  x  0.1 3/(20) 1/2  =  1.13  kcal/mol  from  integrating  over  the 
38.8  A  distance.  Similarly  the  error  of  the  barrier  from  the  DNP 
toward  the  TTF  is  25.5  x  0.13/(14)1/2  =  0.89  kcal/mol  from 
the  integration  over  the  25.5  A  distance.  Hence,  the  small  errors 
in  the  MF  values  can  lead  to  substantial  errors  in  the  PMF  value. 
However,  previous  studies  that  carefully  compare  various  PMF 
calculation  methods  show  that  constraint-biased  sampling  to 
determine  mean  forces  is  one  of  the  best  methods  to  obtain 
reasonable  PMF  values,  even  though,  statistically,  they  contain 
large  error  bars.52 

3.  Results  and  Discussion 

3.1.  Charge  Scheme:  Adiabatic  Approximation.  We  ex¬ 
pected  that  no  set  of  fixed  charges  scheme  would  be  adequate 
enough  to  describe  the  electrostatic  interactions  as  the  highly 
charged  ring  is  moved  along  the  backbone.  Thus,  as  described 
in  section  2.1,  we  obtained  atomic  charges  from  independent 
QM  at  each  position  as  the  ring  is  moved  along  the  backbone. 
This  assumption  of  adiabatically  adjusted  charges  assumes  that 
charge  redistribution  is  much  faster  than  the  time  for  the  ring 
to  travel  along  the  backbone.  To  test  the  effect  of  these  charge 
readjustments  on  the  PMF,  Figure  6  shows  the  PMF  based  on 
three  different  charge  schemes  for  the  neutral  rotaxane  system: 
the  green  curve  was  obtained  using  the  fixed  charges  from  the 
ring  @  TTF,  the  red  curve  was  obtained  using  the  fixed  charges 
from  the  ring  @  DNP,  and  the  blue  curve  was  obtained  using 
adiabatic  charges. 

Clearly,  the  green  and  red  curves  are  biased  to  have  a 
minimum  PMF  at  the  position  for  which  the  charge  was 
calculated,  leading  to  very  bad  estimates  of  the  barrier.  In 
contrast,  the  energy  barrier  between  the  TTF  and  the  DNP  sites, 
based  on  the  adiabatic  charges,  is  consistent  with  experimental 
observations.  Thus,  we  used  the  adiabatic  charges  for  all 
oxidation  states  from  the  neutral  state  to  the  +2  state. 

3.2.  Free  Energy  Profiles  from  PMF  Calculations.  Sam¬ 
pling  the  MFs  from  the  constrained  MD  simulations  (Figure 
7a)  and  integrating  them  along  the  ring  position,  we  calculated 
the  profile  of  the  PMF  for  the  shuttling  motion  of  the  CBPQT4+ 
ring  (Figure  7b).  We  found  that  each  oxidation  state  (neutral 
state  (0),  oxidized  states  (+1  and  +2)),  leads  to  significantly 
different  profiles. 

3.2.1.  AG  T2D-  We  calculated  that  the  most  stable  complex 
for  the  neutral  state  (black)  is  CBPQT4+@TTF  (ring  at  8.9  A) 
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Ring  position  (Angstrom) 

Figure  7.  (a)  Change  of  the  mean  force  as  a  function  of  ring  position 
along  the  backbone,  (b)  Change  of  the  potential  of  mean  force  as  a 
function  of  ring  position  along  the  backbone.  The  green  vertical  lines 
denote  the  ring@TTF  (8.90  A)  and  ring@DNP  (47.70  A). 

whereas  the  CBPQT4+@DNP  state  (ring  at  47.7  A)  is  less  stable 
than  the  TTF  by  AGT2D  =1.0  kcal/mol.  This  agrees  with  various 
experiments,  which  lead  to  AGT2d  —  1.4— 1.6  kcal/mol15,65-68 
on  the  basis  of  the  difference  in  the  binding  free  energies  of 
the  individual  components  of  the  rotaxane  in  the  MeCN  solvent. 
In  addition,  this  calculation  agrees  with  our  previous  computa¬ 
tions  from  QM  (AGT2d  —  2.0  kcal/mol)43  and  Hessian-based 
FF  calculations  (AGT2d  —  2.3  kcal/mol).43  We  emphasize  here 
that  all  previous  experimental  and  theoretical  studies  studied 
AGT2d  by  comparing  the  binding  free  energy  of  the  TTF 
derivatives  with  the  CBPQT4+  ring  and  the  binding  free  energy 
of  the  DNP  derivatives  with  the  CBPQT4+  ring.  Thus,  our 
current  calculation  is  the  first  direct  measurement  of  the  AGT2d 
on  a  complete  rotaxane. 

3.2.2 .  AG  t2d  and  AG  :iD2i  for  Neutral  Case.  We  calculate 
that  the  free  energy  barrier  is  AG*T2d  —  19.03  kcal/mol  from 
the  TTF  toward  the  DNP,  and  AG*D2t  —  18.03  kcal/mol  in  the 
opposite  direction.  The  relaxation  barrier  in  the  neutral  state 
was  measured  for  a  similar  bistable  [2]rotaxane  in  which  the 
triphenylene  spacer  was  missing,  leading  to15  (see  Table  1) 
•AG*e>2t  —  16.2  kcal/mol  (r  ~  0.095  s)  in  MeCN  solvent 
•AG*d2t  —  18.1  kcal/mol  (r  ~  2.4  s)  in  a  MeCN/polymeth- 
ylmethacrylate/propylene  carbonate/LiC104  polymer  matrix 
containing  weight  ratios  of  70:7:20:3 

•AG^£)2t  —  22.21  kcal/mol  (r  ~  2.5  x  103  s)  in  the  molecular- 
switch  tunnel  junction 


In  addition,  the  devices  fabricated  with  this  derivative 
containing  the  triphenylene  spacer  exhibit  a  relaxation  half-life 
of  r  ~  90  min.2  (AG^t  —  22.66  kcal/mol). 

In  addition,  our  free  energy  barrier  is  quite  comparable  to 
the  barriers  to  circumrotation  of  [2]catenanes.  Leigh  and  co¬ 
workers  used  NMR  to  determine  A G*  of  interlocked  catenane 
molecules  as  11—20  kcal/mol  for  various  solvents  and  calculated 
the  free  energy  barrier  as  10—20  kcal/mol  using  force-field- 
based  Hessians.69-71 

Although  our  simulations  were  performed  in  the  gas  phase, 
the  AGW  of  18  kcal/mol  agrees  well  with  the  experimental 
barriers  (17—22  kcal/mol)15,34,49  for  a  variety  of  environments. 
This  suggests  that  the  energy  barrier  does  not  depend  strongly 
on  environment. 

We  did  not  include  the  counterions  in  this  study  because 
preliminary  calculations  showed  that  the  charges  would  some¬ 
times  change  in  erratic  ways  due  to  the  floppy  energy  landscape 
for  the  countercharges.  Indeed,  the  good  agreement  with 
experiment  for  the  barriers  suggests  that  the  instantaneous 
changes  in  the  potential  due  to  counterions  can  be  neglected. 

3.2.3.  AG*tzd  and  AG'  /;?/  for  Oxidized  Cases.  Although  the 
neutral  state  prefers  to  have  the  CBPQT4+  ring  at  the  TTF,  we 
find  that  the  +1  and  +2  oxidized  states  lead  to  a  completely 
different  energy  profile  (Figure  7b).  In  both  cases,  the  DNP 
becomes  the  global  minimum  with  the  TTF  destabilized  by  AG 
=  25.75  kcal/mol  for  the  +1  oxidation  state  and  AG  =  47.78 
kcal/mol  for  the  +2  oxidation  state. 

Starting  with  the  ring  at  the  TTF  site  and  oxidizing,  we  find 
that  the  ring  moves  first  by  ~5  A  to  a  local  minimum  on  the 
ethylene  oxide  linker  (with  an  energy  decrease  by  AG  =  3.25 
kcal  for  the  +1  and  AG  =  5.49  kcal/mol  for  the  +2).  Then,  it 
has  a  free  energy  barrier  of  AG  =  8.70  kcal/mol  (+1  state)  or 
8.02  kcal/mol  (+2  state)  to  continue  past  the  triphenylene  spacer 
and  toward  the  DNP  for  oxidation  states. 

Using  the  Eyring  rate  equation  [Mr  =  (kBT/h)  exp(— AG*/ 
RT)],  the  time  required  to  overcome  this  barrier  to  move  onto 
the  DNP  is  2.9  x  10-7  s  for  the  +1  oxidation  state  and  9.0  x 
10-8  s  for  the  +2  oxidation  state.  It  would  be  interesting  to 
design  an  experiment  to  probe  for  this  predicted  barrier.  It  has 
been  assumed  that  the  huge  Coulomb  potential  of  the  +4  ring 
with  the  +2  TTF  would  preclude  a  barrier.  The  origin  of  this 
barrier  in  the  oxidized  state  is  discussed  below,  which  we  find 
arises  from  the  triphenylene  spacer.  We  expect  that  there  would 
be  no  barrier  without  this  spacer. 

Relative  to  the  final  state  of  the  ring  at  the  DNP  site,  the 
energy  at  the  ethylene  oxide  linker  (EO)  near  the  TTF  site  is 
22.52  kcal/mol  higher  (+1  oxidation),  leading  to  a  Boltzmann 
population  of  10-17.  For  the  +2  oxidation  state,  the  energy  is 
42.41  kcal/mol  higher,  leading  to  a  population  of  10-32.  Thus, 
for  oxidation  states  + 1  and  +2,  we  expect  the  CBPQT4+  ring 
to  stay  on  the  DNP  site  until  the  system  is  reduced. 

Indeed,  there  is  an  experimental  estimate  of  this  reverse  barrier. 
Using  a  modified  AFM  with  the  ring  attached,  Brough  et  al.72 
measured  the  force  exerted  on  the  ring  shuttling  from  the  DNP  to 
the  TTF  in  the  +2  oxidized  system  as  145  pN.  Combining  this 
experimental  data  with  results  from  molecular  mechanics  simula¬ 
tions,  they  estimated  the  energy  barrier  to  be  65  kcal/mol.  This 
can  be  compared  to  our  calculated  barrier  of  50.4  kcal/mol  energy, 
validating  the  accuracy  of  the  experiment.  The  maximum  force 
measured  in  our  simulation  during  the  ring  shuttling  is  583  pN, 
which  is  similar  to  the  experimental  value  of  145  pN. 

3.3.  Effect  of  Coulombic  Energy  and  van  der  Waals 
Energy.  To  understand  why  the  PMF  profiles  are  so  different 
between  the  neutral,  +1,  and  +2  oxidation  states,  we  calculated 
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TABLE  1:  Free  Energy  Barriers,  Rate  Constants,  and  Relaxation  Half-Lives  from  DNP  toward  TTF  (DNP  — *  TTF)  at  298  K 
(All  Simulation  Results  from  this  Work) 

condition 

AG  (kcal/mol)  k  ( s~l)a  zm  (s)“ 


18.03  ±1.5  (simulation) 

0.33  ±  0.83 

2.1  ±5.4 

neutral 

gas  phase  solution 

16.2  ±  0.3  (exp34) 

7.3  ±  3.7 

0.095  ±  0.048 

neutral 

(CH3CN) 

18.1  ±  0.2  (exp15’34’71) 

0.3  ±0.10 

2.4  ±  0.082 

neutral 

polymer  matrix^7 

22.21  ±  0.04  (exp15’34’71) 

(2.7  ±0.19)  x  10“4 

(2.5  ±0.18)  x  103 

neutral 

molecular- switch  junction 

31.22  (simulation) 

6.3  x  10~n 

1.1  x  1010 

oxidation  ±1 

gas  phase 

50.43  (simulation) 

65  kcal/mol  (exp  ±  simulation72) 

4.5  x  10~25 

1.5  x  1024 

oxidation  ±2 
oxidation  ±2 

gas  phase 

SAM  on  Si02  wafers0 

“Values  are  calculated  using  the  Eyring  equation,  1/r  =  ( kBT/h )  exp(— A GXIRT).  ^Weight  ratio 
methacrylate)/propylene  carbonate/LiC104.  c  The  modified  AFM  tip  is  attached  to  the  CBPQT4±  ring. 

70:7:20:3  for  CH3CN/poly(methyl 

Figure  8.  Change  of  Coulombic  interaction  energy  as  a  function  of 
the  ring  position:  (a)  neutral  state;  (b)  oxidation  state  +1;  (c)  oxidation 
state  +2. 


the  change  in  the  Coulombic  interaction  energy  and  the  van 
der  Waals  (vdW)  interaction  energy  as  a  function  of  ring  position 
along  the  backbone,  for  these  three  oxidation  states  (Figures  8 
and  9). 

For  the  neutral  state,  we  find  that  the  Coulombic  energy 
increases  by  60  kcal/mol  as  the  ring  moves  from  the  TTF  to 
the  triphenylene  spacer  (barrier)  and  then  drops  by  45  kcal/mol 
as  it  moves  to  the  DNP.  On  the  other  hand,  the  vdW  energy 
changes,  within  a  range  of  ±4  kcal/mol,  while  the  ring  travels 
from  the  TTF  to  the  DNP. 


Ring  position  (Angstrom) 

Figure  9.  Change  of  van  der  Waals  interaction  energy  as  a  function 
of  the  ring  position:  (a)  neutral  state;  (b)  oxidation  state  +1;  (c) 
oxidation  state  +2. 


Figure  10.  Variations  in  the  total  charge  on  the  ring  as  a  function  of 
the  ring  position  for  the  neutral  case. 


This  indicates  that  the  barrier  is  dominated  by  the  differential 
Coulombic  interactions  with  a  peak  of  443  kcal/mol  at  z  —  28A 
(over  the  spacer).  We  were  quite  surprised  because  we  expected 
the  barrier  to  be  dominated  by  vdW  repulsions  due  to  the  bulky 
size  of  the  triphenylene.  To  understand  why  Coulombic  interac¬ 
tions  are  so  important,  we  plot  in  Figure  10  the  total  charge  on 
the  ring  along  the  pathway,  in  the  neutral  case.  We  see  that  at 
the  TTF  or  the  DNP  positions  there  is  strong  delocalization  from 
the  ring  onto  the  backbone,  but  as  the  ring  passes  over  the 
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triphenylene  spacer  (at  z  —  28  A),  this  charge  localizes  back 
onto  the  ring.  Thus,  we  conclude  that  localization  of  the  ring 
charge  increases  the  Coulombic  repulsion  and  dominates  the 
free  energy  barrier.  This  suggests  that  the  barrier  can  be  modified 
dramatically  by  changing  the  polarity  of  the  spacer. 

We  also  found  that  as  the  system  is  oxidized,  the  magnitude 
of  Coulombic  repulsion  increases  from  380—445  kcal/mol  for 
the  neutral  state,  to  520—545  kcal/mol  for  the  +1  oxidation 
state,  and  finally,  to  670—700  kcal/mol  for  the  +2  oxidation 
state.  In  contrast,  the  vdW  energy  changes  from  165—173  to 
165— 172  to  164—171  kcal/mol  as  the  system  is  oxidized.  This 
implies  that  the  driving  force  inducing  the  mechanical  movement 
of  the  ring  is  the  increased  Coulombic  repulsion  due  to 
oxidization  of  the  rotaxane.  This  confirms  our  view  since  the 
beginning  of  our  experiments. 

However,  the  PMF  profile  (Figure  6b)  still  differs  substan¬ 
tially  from  the  Coulombic  energy  profile  (Figure  8).  For  instance, 
in  the  neutral  case,  the  Coulombic  energy  difference  between 
two  stations  is  12.03  kcal/mol,  which  is  ~12  times  larger  than 
AGT2d  —  1.0  kcal/mol,  and  the  Coulombic  energy  barrier  for 
the  shuttling  from  the  DNP  to  the  TTF  is  46.75  kcal/mol,  which 
is  ~2.6  times  larger  than  AG*D2t  —  18.03  kcal/mol.  Thus,  the 
key  features  of  the  PMF  profile  are  not  fully  explained  in  terms 
of  the  Coulombic  energy  alone.  Another  possible  contributor 
to  the  free  energy  is  vibrational  entropy,  which  can  be 
investigated  directly  from  the  MD  simulation  trajectory.73,74 

4.  Summary 

We  used  constrained  MD  simulations  to  calculate  the  free 
energy  profile  at  300  K  for  the  shuttling  of  the  CBPQT4+  ring 
between  the  TTF  and  the  DNP  in  the  rotaxane  molecule.  This 
free  energy  profile  was  derived  by  calculating  and  integrating 
the  MF  acting  on  the  ring  as  it  is  moved  from  one  position  to 
another  position  along  the  backbone.  We  found  that  it  is 
particularly  important  to  allow  the  charges  to  adjust  adiabatically 
as  the  ring  moves.  Indeed,  we  find  that  the  Coulomb  interactions 
dominate  the  barriers  for  these  systems. 

We  found  that  the  free  energy  barrier  from  the  DNP  to  the 
TTF  is  18.03  kcal/mol  for  the  neutral  system,  which  agrees  well 
with  experimental  values  of  17—22  kcal/mol  for  various 
environments.  We  calculate  that  the  AG  between  the  TTF  and 
the  DNP  positions  is  1.0  kcal/mol,  which  compares  well  with 
experimental  results  of  1.4— 1.6  kcal/mol  obtained  from  binding 
energies  of  separate  DNP  and  TTF  systems  with  the  CBPQT4+ 
ring. 

These  results  validate  the  accuracy  of  our  computational 
procedure.  Thus,  we  can  now  use  this  validated  technique  for 
estimating  the  switching  kinetics  for  new  designs  of  molecular 
architectures. 
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